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Abstract. We study the spin dynamics of a long nanomagnet driven by an 
electrical current. In the case of only DC current, the spin dynamics has a 
sophisticated bifurcation diagram of attractors. One type of attractors is a 
weak chaos. On the other hand, in the case of only AC current, the spin 
dynamics has a rather simple bifurcation diagram of attractors. That is, for 
small Gilbert damping, when the AC current is below a critical value, the 
attractor is a limit cycle; above the critical value, the attractor is chaotic 
(turbulent). For normal Gilbert damping, the attractor is always a limit cycle 
in the physically interesting range of the AC current. We also developed 
a Melnikov integral theory for a theoretical prediction on the occurrence of 
chaos. Our Melnikov prediction seems performing quite well in the DC case. 
In the AC case, our Melnikov prediction seems predicting transient chaos. 
The sustained chaotic attractor seems to have extra support from parametric 
resonance leading to a turbulent state. 
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1. Introduction 

The greatest potential of the theory of chaos in partial differential equations 
lies in its abundant applications in science and engineering. The variety of the spe- 
cific problems demands continuing innovation of the theory [17] [16] [18] [23] [24] 
[15j [20] [21j . In these representative publications, two theories were developed. 
The theory developed in [17j |16j |18] involves transversal homoclinic orbits, and 
shadowing technique is used to prove the existence of chaos. This theory is very 
complete. The theory in |23] |24) [15] [20] [21] deals with Silnikov homoclinic 
orbits, and geometric construction of Smale horseshoes is employed. This theory 
is not very complete. The main machineries for locating homoclinic orbits are (1). 
Darboux transformations, (2). Isospectral theory, (3). Persistence of invariant man- 
ifolds and Fenichel fibers, (4). Melnikov analysis and shooting technique. Overall, 
the two theories on chaos in partial differential equations are results of combining 
Integrable Theory, Dynamical System Theory, and Partial Differential Equations 

In this article, we are interested in the chaotic spin dynamics in a long nano- 
magnet diven by an electrical current. We hope that the abundant spin dynamics 
revealed by this study can generate experimental studies on long nanomagets. To 
illustrate the general significance of the spin dynamics, in particular the magneti- 
zation reversal issue, we use a daily example: The memory of the hard drive of a 
computer. The magnetization is polarized along the direction of the external mag- 
netic field. By reversing the external magnetic field, magnetization reversal can 
be accomplished; thereby, generating and 1 binary sequence and accomplishing 
memory purpose. Memory capacity and speed via such a technique have reached 
their limits. The "bit" writing scheme based on such Oersted-Maxwell magnetic 
field (generated by an electrical current) encounters fundamental problem from clas- 
sical electromagnetism: the long range magnetic field leads to unwanted writing or 
erasing of closely packed neighboring magnetic elements in the extremely high den- 
sity memory device and the induction laws place an upper limit on the memory 
speed due to slow rise-and-decay-time imposed by the law of induction. Discovered 
by Slonczewski |35j and Berger [l], electrical current can directly apply a large 
torque to a ferromagnet. If electrical current can be directly applied to achieve 
magnetization reversal, such a technique will dramatically increase the memory 
capacity and speed of a hard drive. The magnetization can then be switched 
on the scale of nanoseconds and nanometers [39j . The industrial value will be 
tremendous. Nanomagnets driven by currents has been intensively studied recently 

[laiioiiMiiiiiiaiiiiiMiiiiiMiiiiiiaaisaiMiiiia The 

researches have gone beyond the original spin valve system |35j [l]. For instance, 
current driven torques have been applied to magnetic tunnel junctions |36j [6], di- 
lute magnetic semiconductors [40] . multi-magnet couplings |10) |14] . AC currents 
were also applied to generate spin torque [34, [7^ . Such AC current can be used to 
generate the external magnetic field [34] or applied directly to generate spin torque 

Mathematically, the electrical current introduces a spin torque forcing term 
in the conventional Landau-Lifshitz-Gilbert (LLC) equation. The AC current can 
induce novel dynamics of the LLC equation, like synchronization [34] [7] and chaos 
[25] |41| . Both synchronization and chaos are important phenomena to understand 
before implementing the memory technology. In [25] |41j , we studied the dynamics 
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of synchronization and chaos for the LLG equation by ignoring the exchange field 
(i.e. LLG ordinary diH'erential equations). When the nanomagnetic device has the 
same order of length along every direction, exchange field is not important, and we 
have a so-called single domain situation where the spin dynamics is governed by 
the LLG ordinary differential equations. In this article, we will study what we call 
"long nanomagnet" which is much longer along one direction than other directions. 
In such a situation, the exchange field will be important. This leads to a LLG 
partial differential equation. In fact, we will study the case where the exchange 
field plays a dominant role. 

The article is organized as follows: Section 2 presents the mathematical formu- 
lation of the problem. Section 3 is an integrable study on the Heisenberg equation. 
Based upon Section 3, Section 4 builds the Melnikov integral theory for predicting 
chaos. Section 5 presents the numerical simulations. Section 6 is an appendix to 
Section 3. 

2. Mathematical Formulation 

To simplify the study, we will investigate the case that the magnetization de- 
pends on only one spatial variable, and has periodic boundary condition in this 
spatial variable. The application of this situation will be a large ring shape nano- 
magnet. Thus, we shall study the following forced Landau-Lifshitz-Gilbert (LLG) 
equation in the dimensionless form, 

(2.1) dtm = -m x H - earn x {m x H) + e(/3i -|- P2 cosojot)m x (m x ex) , 
subject to the periodic boundary condition 

(2.2) m{t, x + 2'k)= m{t, x) , 

where m is a unit magnetization vector m, = (mi, TO2, TO3) in which the three 
components are along (x, y, z) directions with unit vectors (e^, e^, e^), \m\{t, x) = 1, 
the effective magnetic field H has several terms 

^ = ^exch + ^ext + -f^dem + -'^ani 

(2.3) = d^m + eaex — em^ez + ebmiex , 

where ^fgxch ~ ^s"^ exchange field, H^^-^ = eae^ is the external field, 

■'^dem ~ — ewsCz is the demagnetization field, and -f^a,ni ~ ebmiex is the anisotropy 
field. For the materials of the experimental interest, the dimensionless parameters 
are in the ranges 

a « 0.05, 6 « 0.025, a « 0.02 , 

(2.4) /3i e [0.01,0.3] , /32 G [0.01,0.3] ; 

and e is a small parameter measuring the length scale of the exchange field. One 
can also add an AC current effect in the external field -ffe^t , but the results on the 
dynamics are similar. 

Our goal is to build a Melnikov function for the LLG equation around domain 
walls. The roots of such a Melnikov function provide a good indication of chaos. 

For the rest of this section, we will introduce a few interesting notations. The 
Pauli matrices are: 
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Let 

(2.6) m+ ~ mi + im2 , to_ = mi — im2 , 

i.e. fn^ ~ m_. Let 



Hs Hi - 1H2 

Hi + tH2 -H3 



(2.7) 1 — miCTi = 

J J \ rn^ —m^ 

Thus, = / (the identity matrix). Let 

H = —H — am x H + (3m x Cx , H = 
Then the LLG can be written in the form 

(2.8) z9tr = i[r,n], 

where [r,n] = FH - HL. 

3. Isospectral Integrable Theory for the Heisenberg Equation 

Setting e to zero, the LLG (|2.ip reduces to the Heisenberg ferromagnet equation, 

(3.1) dtm — —711 X m^x ■ 

Using the matrix F introduced in (|2.7p , the Heisenberg equation p.ip has the form 

(3.2) idtV^-Uv^Yxx], 



where the bracket [ , ] is defined in ()2.8p . Obvious constants of motion of the 
Heisenberg equation (j3.ip are the Hamiltonian, 



2tt 



dx 







the momentum, 



and the total spin. 



mim2x - m2mix , 

ax 

1 + ma 



27r 

mdx . 







The Heisenberg equation p.ip is an integrable system with the following Lax pair, 

(3.3) dxip = iXTip , 

(3.4) dtij = -^(4zAr+[F,r,])V , 

where ip = {ipi, 1^2)'^ is complex- valued, A is a complex parameter, F is the matrix 
defined in (|2.7p . and [F,Fj;] ~ TTx — TxT. In fact, there is a connection be- 
tween the Heisenberg equation p.ip and the ID integrable focusing cubic nonlinear 
Schrodinger (NLS) equation via a nontrivial gauge transformation. The details of 
this connection are given in the Appendix. 
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3.1. A Simple Linear Stability Calculation. As shown in the Appendix, 
the temporaUy periodic solutions of the NLS equation correspond to the domain 
wahs of the Heisenberg equation. Consider the domain wall 

Fq = ( ^j^j, Q 1 , ^ e Z; i.e. mi = cos^x, m2 = sin^x, ma = 0, 

which is a fixed point of the Heisenberg equation. Linearizing the Heisenberg equa- 
tion at this fixed point, one gets 



Let 



we get 



Let 



r = 



ma e *^^(mi — ?m2) 

e'^^(mi + im2) —ma 



dtuii = , dtm2 = nisxx + S, ms , dtm^ = -m2xx - 2$mia; 



'J — 

k=0 



(m^^, (t) COS kx + m-f, (t) sin kx) 



where 'rn%{t) = c^^e^*, c% and fl are constants. We obtain that 



(3.5) fl^ y/k^{^^-k^) 

which shows that only the modes < < |^| are unstable. Such instability is 
called a modulational instability, also called a side-band instability. Comparing 
the Heisenberg ferromagnet equation ()3.ip and the Landau-Lifshitz-Gilbert equa- 
tion ()2.ip . we see that if we drop the exchange field -ffg^ch ~ ^^"^ the effective 
magnetic field H (|2.3|) . such a modulational instability will disappear, and the 
Landau-Lifshitz-Gilbert equation (|2.f p reduces to a system of three ordinary differ- 
ential equations, which has no chaos as verified numerically. Thus the modulational 
instability is the source of the chaotic magnetization dynamics. 
In terms of rri^f., we have 

J^^ik = 0. J^^2k = (C' - k^)^ik^ J^^tk = fc'"*^fc T 2^fcmf^. 
Choosing ^ = 2, we have for fc = 0, 



for fc = 1, 





Cl 







■:•) 


+ C2 1 V3 








for k 










+ C3 1 


T8i / 


V 4i 
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for fc> 2, 

1 \ / 

ci ( ±4/fc +C2 VA:2-4cos(fcVA:2-4t) 
/ \ fcsin(fcVfc2 - 4i) 



+C3 ( -Vfc^ - 4sin(fc\/A;2 - At) 
kcos{ky/k^ - At) 

where ci, C2 and C3 are arbitrary constants. 

The nonhnear foUation of the above hnear modulational instabiUty can be es- 
tabUshed via a Darboux transformation. 

3.2. A Darboux Transformation. A Darboux transformation for (|3.3p - 
can be obtained. 



Theorem 3.1. Let (j) = (0i,02)^ be a solution to the Lax pair l^.g)J - P^ at 
(r,!/). Define the matrix 

where 

^ h ~4>2 



N 

' "2 <Pi 

Then if solves the Lax pair i3.3\ )- (K4^ at (T, A^, 
(3.6) iA = G^j 

solves the Lax pair i3.3\) - (K4\ l at (T, X), where T is given by 



(3.7) r = iV ^ ^ J N-^TN ^ ^ ^„., J iV 

where e*^ = 

The transformation (|3.6p - (|3.7p is called a Darboux transformation. This theo- 
rem can be proved either through the connection between the Heisenberg equation 
and the NLS equation (with a well-known Darboux transformation) f2| , or through 
a direct calculation. 

Notice also that = /. Let 



$1 -$2 A f e~ 



-1 



(3.8) 
Then 



$2 t)-''[0 e-)^ 
|<^iP + |'^2P V -2ismeM2 e''^|0i|2-He"^«|,^2p 



(3.9) r = 

where 



ma mi — im2 

mi -|- im2 —rn-s 



m_)_ = mi + im2 = $1 (mi + im2) — $2('^i ~ ^"^2) + 2<i>i$2m'3 , 
= d^"!!^ - 1*2!^) m3 - $i$2(mi -I- im2) - $i$2(™i - im2) ■ 
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One can generate the figure eight connecting to the domain waU, as the nonhnear 
fohation of the modulational instabihty, via the above Darboux transfomation. 



3.3. Figure Eight Connecting to the Domain Wall . Let F be the domain 



wah 



r = 







-i2x 







i.e. mi = cos 2a;, m2 = sin 2a;, and 777,3 = 0. Solving the Lax pair (|3.3p - (|3.4p . one 
gets two Bloch eigenfunctions 

2Aexp{i(A;-2)a;} 
(fc-2)exp{i(fc + 2)a;} 

To apply the Darboux transformation (j3.7[) , we start with the two Bloch functions 
with k — ±1, 



(3.10) ^ 



nt 



-iXk 



k = ±2^1 + A2 . 



(3.11) 



The wise choice for 



\/3e- 



—le ' 
V3e*' 



exp 



exp ■ 




V3, 



(3.12) 



used in (|3.7p is: 



where c+/c = exp{fT + 77}, r = ^{V^t + a), and x - 
Darboux transformation (|3.7p . one gets 

2 sech2rcos2x 



ie^+'x + V3e-^-*x) e 
1 



[X 



7). Then from the 



7771 + 77772 



(3.13) 

(3.14) 
As t ^ 



7773 



(2 - V3 sech2Tsin2x)2 

+7 (V3-2 sech2r sin 2x 

2 sech2T tanh 2t cos 2x 
(2- V3sech2rsin2x)2 



sech2r cos 2x 



±00, 



TTT-i 







- cos 2a; , 7772 — s- — sin 2x , 7773 
The expressions (|3.13l) - (|3.14l) represent the two dimensional figure eight separatrix 
connecting to the domain wall (777+ = — e'^^, 7773 = 0), parametrized by a and 
7. See Figure [1] for an illustration. Choosing 7 = 0, tt, one gets the figure eight 
curve section of Figure [T] The spatial-temporal profiles corresponding to the two 
lobes of the figure eight curve are shown in Figure [2l In fact, the two profiles 
corresponding the two lobes are spatial translates of each other by tt. Inside one 
of the lobe, the spatial-temporal profile is shown in Figure [31[a). Outside the figure 
eight curve, the spatial-temporal profile is shown in Figure [3jb). Here the inside 
and outside spatial-temporal profiles are calculated by using the integrable finite 
difference discretization [9] of the Heisenberg equation (j3.1|) . 



(3.15) 



dt 



m{j) 



m{j) X 



m{j + 1) 



Mj - 1) 



1 -I- 777(j) • 777(j -|- 1) 1 + 777(j — 1) • m{j) ^ 
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where m(j) = m{t,jh), j = 1, - ■ ■ ,N, Nh — 2tt, and h is the spatial mesh size. For 
the computation of Figm'e[3l we choose N — 128. 
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Figure 1. The separatrix connecting to the domain waU m+ = 




(a) 7 = (b) 7 = TT 



Figure 2 . The spatial-temporal profiles corresponding to the two 
lobes of the figure eight curve. 

By a translation x ^ x + 9, one can generate a circle of domain walls: 

m+ = -e*2(^+«) , mg = , 

where 9 is the phase parameter. The three dimensional figure eight separatrix 
connecting to the circle of domain walls, parametrized by ct, 7 and 9] is illustrated 
in Figure |4l 

In general, the unimodal equilibrium manifold can be sought as follows: Let 

nij — Cj cos 2x + Sj sin 2x , j — 1, 2, 3, 

then the uni-length condition |m|(x) = 1 leads to 

|c| = 1 , |s| = 1 , c - s = , 

where c and s are the two vectors with components Cj and Sj . Thus the unimodal 
equilibrium manifold is three dimensional and can be represented as in Figure \5\ 
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(a) inside (b) outside 



Figure 3. The spatial-temporal profiles corresponding to the in- 
side and outside of the figure eight curve. 




Figure 4. The separatrix connecting to the circle of domain walls 

Using the formulae (|3.13[) - (|3.14[) . we want to build a Melnikov integral. The 
zeros of the Melnikov integral will give a prediction on the existence of chaos. To 
build such a Melnikov integral, we need to first develop a Melnikov vector. This 
requires Floquet theory of l|3.3p . 

3.4. Floquet Theory. Focusing on the spatial part (|3.3p of the Lax pair 
let Y{x) be the fundamental matrix solution of Y{0) = / (2 x 2 

identity matrix), then the Floquet discriminant is defined by 

A = trace y(27r) . 

The Floquet spectrum is given by 

cr = {XeC\ -2< A(A) < 2} . 
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Figure 5. A representation of the 3 dimensional unimodal equi- 
librium manifold. 

Periodic and anti-periodic points (which correspond to periodic and anti-periodic 
eigenfunctions respectively) are defined by 

A(A±) = ±2 . 

A critical point A^'^^ is defined by 

A multiple point A*-"-* is a periodic or anti-periodic point which is also a critical 
point. The algebraic multiplicity of A^"^ is defined as the order of the zero of 
A(A) ± 2 at A("). When the order is 2, we call the multiple point a double point, 
and denote it by A*^''-'. The order can exceed 2. The geometric multiplicity of A*^"^ 
is defined as the dimension of the periodic or anti-periodic eigenspace at A^"-*, and 
is either 1 or 2. 

Counting lemmas for A* and A'^'^^ can be established as in [30] [23] . which lead 
to the existence of the sequences {A^} and {Aj*^^} and their approximate locations. 
Nevertheless, counting lemmas are not necessary here. For any A G C, A(A) is a 
constant of motion of the Heisenberg equation p.ip . This is the so-called isospectral 
theory. 

Example 3.2. For the domain wall toi — cos2x, m2 = sin 2a;, and ma — 0; the 
two Bloch eigenfunctions are given in (j3.10p . The Floquet discriminant is given by 



A = 2 cos 
The periodic points are given by 

A 



27rVl + A2 



± 



1 



n G Z , 



n IS even 



The anti-periodic points are given by 
A 



± 



1 , n G Z , n is odd 



The choice of g!)+ and 



A' 



correspond to n = ±1 and A 
A 



jV3/2 with fc = ±1. 



-47r- 



Vl + A^ 



A" = -47r(l-f A2)-3/2sin 27rVT+A2 



27rVl + A2 

o 2 A' 

-Svr^- — 

1 A2 



27rVl + A2 
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lAj 



-e — e- 



-e — e- 



FiGURE 6. The periodic and anti-periodic points corresponding to 
the potential of domain wall m+ = e*^^, = 0. The open circles 
are double points, the solid circle at the origin is a multiple point 
of order 4, and the two bars intersect the imaginary axis at two 
periodic points which are not critical points. 



When n = 0, i.e. Vl + = 0, by L'Hospital's rule 

A' ^ -Svr^A , X^±i . 

That is, X = ±i are periodic points, not critical points. When n — ±1, we have two 
imaginary double points 

A = ±iV3/2 . 

When n = ±2, A = is a multiple point of order 4. The rest periodic and anti- 
periodic points are all real double points. FigurelHlis an illustration of these spectral 
points. 



3.5. Melnikov Vectors. Starting from the Floquet theory, one can build 
Melnikov vectors. 

Definition 3.3. An important sequence of invariants Fj of the Heisenberg 
equation is defined by 



Lemma 3.4. If{X) '} is a simple critical point of A, then 

dFj _ dA 
dm dm 



Proof. We know that 

dF^ 
dm dm 



dFj _ dA 



dA 
5A 



dX 



dm 



Since 



we have 



dA 
9A 



= , 



A=A'.' 



d^A 



9A2 



dX 



d^A 



A=A 



■c) dm dXdr. 



= 



a=a(^' 
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(c) 

Since A,- is a simple critical point of A, 



Thus 









^0. 














d^A 


dm 


9A2 




dXdm 



Notice that A is an entire function of A and m |23 , then we know that 



bounded, and 



dm 



dFj _ dA 
dm dm 



IS 



□ 



Theorem 3.5. As a function of two variables, A — A(A, m) has the partial 
derivatives given by Bloch functions tjj^ (i.e. 'il'^{x) = e^-^^^jj^ {x) , where V'^ ttTC 
periodic in x of period 2tt, and A is a complex constant): 

dA V^—A 

dA 
dm- 
dA 
dms 
dA 

~dx ~ V(V'+,V'-)7o 

where W{ip'^ ,1^^) ~ ^ "4^2 "^i '^^ Wronskian. 



VA2 -4 



Proof. Recall that Y is the fundamental matrix solution of (|3.3[) . we have the 
equation for the differential of Y 

d.j:dY ^ iXTdY + i{dXT + Xdr)Y , dY{0) = 0. 

Using the method of variation of parameters, we let 

dY = YQ , Q(0) = . 

Thus 

Q{x)^i I Y-\dXr + XdT)Ydx , 
Jo 

and 



Finally 

(3.16) 

Let 



dY{x) =iY [ Y-\dXT + XdT)Ydx . 
Jo 

dA = trace dY {2n) 

= z trace |r(27r)j^ ^-^(dAr + Adr)rda;| 
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where are two hnearly independent Bloch functions (For the case that there is 
only one hnearly independent Bloch function, L'Hospital's rule has to be used, for 
details, see (23 ), such that 

where are periodic in x of period 27r and A is a complex constant (The existence 
of such functions is the result of the well known Floquct theorem). Then 

Z{x) = Y{x)Z{0) , Y{x) = Z{x)[Z{0)]-^ . 



Notice that 

Then 
Thus 
and 



Z{2tt) = Z{0)E , where E = 



e^'^ 

n „-A27r 



y(27r) ^ Z{{))E[Z{Q)]-^ . 
A = trace y(27r) = trace E = e^^^r g-A27r ^ 

g±A2^ ^ i[A± -4] . 
In terms of Z , dA as given in p.l6p takes the form 

dA = i trace |z(0)£;[Z(0)]"^^ Z{Q)[Z{x)]-^ {d\T + XdT)Z{x)[Z{Q)]-^ dx^ 
= i trace (i;/ [Z{xyf^{d\T + \dT)Z{x)dx\ 



from which one obtains the partial derivatives of A as stated in the theorem. □ 

It turns out that the partial derivatives of Fj provide the perfect Melnikov 
vectors rather than those of the Hamiltonian or other invariants |23j , in the sense 
that Fj is the invariant whose level sets are the separatriccs. 

3.6. An Explicit Expression of the Melnikov Vector Along the Figure 
Eight Connecting to the Domain Wall. We continue the calculation in sub- 
section 13.31 to obtain an explicit expression of the Melnikov vector along the figure 
eight connecting to the domain wall. Apply the Darboux transformation (|3.6p to 
(jy^ (|3.1ip at A = t^, we obtain 

7± _ ^ v-vey.^{T \a T i\l}W{(j)^ 



In the formula (|3.6p . for general A, 

detG = -^^^-^44^^-^ , Wi'i+.'ip-) = delGWi'tP+.ip-) . 

In a neighborhood of A = j^, 

A^ - 4 = A(i^)A"(i/)(A - vf + higher order terms in (A - v) 
As A ^ j^, by L'Hospital's rule 



VA2 -4 ^ I\{v)IS!'{y) 
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Notice, by the calculation in Example 13.21 that 



then by Theorem 13.51 
9A 



dA 



A{iy) = -2 , A"(zy) = -247r2 , 



9m_ 
dA 



dra^ 



= 12\/37r 
= 12%/37r 



2i 



where rh is given in p.l3p - p.l4p . With the explicit expression (j3.12l) of 0, we 
obtain the explicit expressions of the Melnikov vector, 



dA 



dm A 



i sech2T 



(3.17) 



dA 



dm- 



2 (2- x/3 sech2rsin2x)2 
+i(2 - tanh 2r) sin 2x - iVs sech2T 
Sa/Stt —i sech2r 



(1 - 2tanh2T) cos2x 



(3.18) 



dA 



dms 

(3.19) 

where again 



2 (2-\/3sech2rsin2x)2 
-i{2 + tanh 2r) sin 2x + iVs sech2T 
3-\/37r 2i sech2T 



(1 + 2tanh2T) cos2x 



i2x 



2 (2-%/3sech2rsin2x)2 
-iVS tanh 2t cos 2x 



2 sech2T — \/3sin2x 



\/3 cr 1 , . 



and (T and 7 are two real parameters. 

4. A Melnikov Function 

The forced Landau-Lifshitz-Gilbcrt (LLG) equation ()2.ip can be rewritten in 
the form, 

(4.1) dtm = -m X m^x + f^f + f^^g 

where / is the perturbation 

/ = —am X Cx + 'm^[m x Cz) — hmiijn x Cx) 

—am X (m x mxx) + {Pi + /?2 cosujQt)m x (m x ex) , 

g ~ —am x [to x {aCx ~ m^Cz + bmiCx)] ■ 
The Melnikov function for the forced LLG (|2.ip is given as 



M 



00 p27T " 



00 JO 



■7^ ifi + 1/2) + 7^ (/i - 4/2) + J. — /s 

OTO+ OTO_ OTOa 



dxdt 
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where m is given in p.l3p - (l3.14|) . and ^ {w — m+, m_, 7713) are given in (|3.17p - 
p.l9p . The Melnikov function depends on several external and internal parameters 
M = M{a, b, a, /32, i^Oj f , 7) where a and 7 are internal parameters. We can split 
/ as follows: 



where 



+P2 



/3i/ 



(/3i 



I 



is) - 



-m X ex , 
m,z{m X e^) , 
-mi{m X e^) , 
— m X (rn x mj;^;) , 
m X (m X Bj;) , 



cos ( --j=u}qt j m X (m x e^;) 



sin ( --j=ujqt j TO X (to X e^;) 



Thus M can be splitted as 



M 



(4.2) 



+f32 



cos I —r^ I M^'') 
V V3, 



sin 

lV3 



■ sm I I 



where M(f) = M(^)(7), C = a,0,6,a,/3i, and Aff^) = Af^^) (7, wq), C = c, s. 

In general [19| . the zeros of the Melnikov function indicate the intersection of 
certain center-unstable and center-stable manifolds. In fact, the Melnikov function 
is the leading order term of the distance between the center-unstable and center- 
stable manifolds. In some cases, such an intersection can lead to homoclinic orbits 
and homoclinic chaos. Here in the current problem, we do not have an invariant 
manifold result. Therefore, our calculation on the Melnikov function is purely from 
a physics, rather than rigorous mathematics, point of view. 

In terms of the variables to+ and 7713, the forced Landau-Lifshitz-Gilbert (LLG) 
equation (|2.ip can be rewritten in the form that will be more convenient for the 
calculation of the Melnikov function, 

(4.3) dtm+ = i{m+m:ixx - m^m+xx) + e/+ + e^5+ . 

(4.4) dt-m^ = ^{m+Trq:^^-rn^m+xx) + <^f3 + <^'^93, 
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/+ = 



/3 = 



9+ 
93 



+P2 



91 + i92 = aag^^^ + ag^°^ + abg^^^ , 



-ims 



= —11713171+ , 

= -i^m3{m++rfq:) , 
= ^rn+{m+-m^)-ml , 



cos I —^uJoT 



sm 



Wot 



m+(m_|_ — m+) — 
m+ (m+ — rrT+ ) — 



/f = 0, 

/^"' = rn3xx|?n+|^-im3 (m+rrEj:^^ + m:pm+ ^ 



/F'^ = ^"^3(m+ +m+) , 
/s''^ = ^ cos ^-^woT^ m3(m+ + m^) , 
/s*"* = ^ sin ^;^a;or^ m3(m+ + rnT) , 
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(a) 

9+ 



1-+) J 



(0) 
9+ 

J'') 



(a) 
53 

(0) 
53 

(b) 
53 



1 , , 

--TO3(m+ + TO+) , 



Direct calculation gives that 

M('')(7) = M(°)(7) = A/('')(7) = 0, M(")(7) = 91.3343, 

and M^'^i) and M^^^ are real, while Af*^^) is imaginary. The graph of Af^'^^^ is 
shown in Figure[7lja) (Notice that Af'^^' is independent of ujq). The graph of A/'-'^^ 
is shown in Figure [71(b). The imaginary part of Af'^'*^ is shown in Figure [Tjc). In 
the case of only DC current (/?2 = 0), Af = (|4?2l) leads to 

(4.5) a = A/(^i) /91.3343, 

where M^^'-'^ {■-■/) is a function of the internal parameter 7 as shown in Figure [Tlja). 
In the general case {P2 ^ 0), Af^*^ (7,1^0) — determines curves 

(4.6) 7 = 7(c^o) = 0, 7r/2, tt, 3^2, 
and M = IUSI leads to 



(4.7) > 



(91.3343a + /?iM('3i) 



cos 



where M^^^i) and Af('=) are evaluated along the curve (gj]), Af('^i) = ±43.858 ('+' 
for 7 = 7r/2, 37r/2; '— ' for 7 = 0, tt), A/^^^ is plotted in Figure [5] (upper curve 
corresponds to 7 = 7r/2, 37r/2; lower curve corresponds to 7 = 0, tt), and 

fauJo\ _ 91.3343a + ^iAf('5i) 

5. Numerical Simulation 

In the entire article, we use the finite difference method to numerically simulate 
the LLG (|2.ip . Due to an integrable discretization [9| of the Heisenberg equation 
p.ip . the finite difference performs much better than Galerkin Fourier mode trun- 
cations. As in (|3.15p . let m(j) = 'm{t,jh), j — 1, ■ ■ ■ , N, Nh — 2tt, and h is the 
spatial mesh size. Without further notice, we always choose N = 128 (which pro- 
vides enough precision). The only tricky part in the finite difference discretization 
of (12. ip is the second derivative term in H, for the rest terms, just evaluate m at 
m{j): 

ct2 f 2 1^ m{j + 1) m{j - 1) 



/i^ \ 1 + m(j) • m(j + 1) 1 + m(j — 1) • m(j) 
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Figure 7. (a). The graph of M^^i) as a function of 7, and M^^^i) 
is independent of ujq- (b). The graph of M^'^' as a function of 7 
and Wo- (c). The graph of the imaginary part of M^") as a function 

of 7 and luq. 

5.1. Only DC Current Case. In this case, /32 = in (|2.ip . and we choose 
(3i as the bifurcation parameter, and the rest parameters as: 

(5.1) a = 0.05, b = 0.025, a = 0.02, e = 0.01. 

The computation is first run for the time interval [0, 81207r], then the figures are 
plotted starting from t = 81207r. The bifurcation diagram for the attractors, and 
the typical spatial profiles on the attractors are shown in Figure [SI This figure 
indiactes that interesting bifurcations happen over the interval /3i £ [0,0.15] which 
is the physically important regime where Pi is comparable with values of other 
parameters. There are six bifurcation thresholds ci • • • cg (Fig. [5]). When /3i < ci, 
the attractor is the spatially uniform fixed point mi — 1 (m2 — = 0). When 
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I I II I III I IV I V I VI I VII 



Figure 9. The bifurcation diagram for the attractors and typ- 
ical spatial profiles on the attractors in the case of only DC 
current where Pi is the bifurcation parameter, ci • • • cg are the 
bifurcation thresholds, and ci G [0,0.01], C2 G [0.0205,0.021], 
C3 G [0.0231,0.0232], C4 G [0.025,0.026], C5 G [0.08,0.1], cg G 
[0.13,0.15]. 



ci < /3i < C2 , the attractor is a spatially non- uniform fixed point as shown in Figure 
llOf aV When C2 < /3i < C3, the attractor is spatially non- uniform and temporally 
periodic (a limit cycle) or quasiperiodic (a limit torus) as shown in Figure fTOTb). 
Here as the value of /3i is increased, first there is one basic temporal frequence, 
then more frequencies enter and the temporal oscillation amplitude becomes bigger 
and bigger. When C3 < /3i < C4, the attractor is chaotic, i.e. a strange attractor 
as shown in Figure llOf c) . Even though the chaotic nature is not very apparent 
in Figure [TOFc). due to the smallness of the perturbation parameter together with 
smallness of all other parameters, the temporal evolution is chaotic, and we have 
used Liapunov exponent and power spectrum devices to verify this. When C4 < 
/3i < C5, the attractor is spatially non-uniform and temporally periodic (a limit 
cycle) as shown in Figure ITIT a). The spatial modulation is small, and it is even not 
apparent in Figure fTlT a). But it is apparent on the individual typical spatial profiles 
as seen in region V in Figure [9l With the increase of f3i , the spatial modulation 
becomes smaller and smaller. When C5 < /3i < cq, the attractor is spatially uniform 
and temporally periodic (a limit cycle, the so-called procession) as shown in Figure 
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(a) /3i = 0.0205 spatially non-uniform fixed (b) /3i = 0.023 spatially non-uniform and 
point temporally periodic or quasiperiodic attrac- 

tor 




(c) /3i = 0.0235 weak chaotic attractor 

Figure 10. The spatio-temporal profiles of solutions in the at- 
tractors in the case of only DC current. 



IllT b). When f3i > cg, the attractor is the spatially uniform fixed point mi = — 1 
(to2 — = 0)- 

When (32 — 0, the Melnikov function predicts that around /3i — 0.041 (|4.5p . 
there is probably chaos; while the numerical calculation shows that there is an 
interval [0.0231,0.026] for f3i where chaos is the attractor. Since the perturbation 
parameter e = 0.01 in the numerical calculation, the Melnikov function prediction 
seems in agreement with the numerical calculation. 

Some of the attractors in Figure [S] are attractors of the corresponding ordinary 
differential equations by setting dx ~ in (|2.ip . i.e. the single domain case. These 
attractors are the ones in regions I, VI and VII in Figure [S] |25j |41j . Because we 
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(a) /3i = 0.026 spatially non-uniform and (b) /3i = 0.1 spatially uniform and tempo- 
tcmporally periodic attractor rally periodic attractor 



Figure 11. The spatio-temporal profiles of solutions in the at- 
tractors in the case of only DC current (continued). 

are studying the only DC current case, the ordinary differential equations do not 
have any chaotic attractor [25] [41j . 

On the other hand, the partial differential equations (|2.1[1 does have a chaotic 
attractor (region IV in Figure (5]). 

In general, when /3i < 0, the Gilbert damping dominates the spin torque driven 
by DC current and mi = 1 is the attractor. When /3i > 0.15, the spin torque driven 
by DC current dominates the Gilbert damping, toi = — 1 is the attractor, and we 
have a magnetization reversal. In some technological applications, /3i > 0.15 may 
correspond too high DC current that can burn the device. On the other hand, in 
the technologically advantageous interval /3i e [0, 0.15], magnetization reversal may 
be hard to achieve due to the sophisticated bifurcations in Figure O 

5.2. Only AC Current Case. In this case, /3i = in (|2.ip . and we choose 
(32 as the bifurcation parameter, and the rest parameters as: 

(5.2) a = 0.05, b = 0.025, a = 0.0015, e = 0.01, ujq = 0.2. 

Unlike the DC case, here the figures are plotted starting from i = 0. It turns 
out that the types of attractors in the AC case are simpler than those of the 
DC case. When (32 = 0, the attractor is a spatially non-uniform fixed point as 
shown in Figure [121 In this case, the only perturbation is the Gilbert damping 
which damps the evolution to such a fixed point. When Q < (32 < P2 where 
(32 G [0.18,0.19], the attractor is a spatially non-uniform and temporally periodic 
solution. When (32 > /SJ ; the attractor is chaotic as shown in Figure [TH Our 
Melnikov prediction (|4.7p predicts that when \(32\ > 0.003, certain center- unstable 
and center-stable manifolds intersect. Our numerics shows that such an intersection 
seems leading to transient chaos. Only when \(32\ > (^2^ the chaos can be sustained 
as an attractor. It seems that such sustained chaotic attractor gains extra support 
from parametric resonance due to the AC current driving |22] , as can be seen from 
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"Ox t 
(a) /32 = spatially non-uniform fixed point (b) /32 = temporal evolution at one spatial 



location 




"Ox t 

(c) /32 = 0.21 chaotic attractor (d) = 0.21 temporal evolution at one spa- 

tial location 



Figure 12. The attractors in the case of only AC current. 



the turbulent spatial structure of the chaotic attractor (Fig. [121), which diverges 
quite far away from the initial condition. Another factor that may be relevant 
is the fact that higher-frequency spatially oscillating domain walls have more and 
stronger linearly unstable modes (|3.5[) . By properly choosing initial conditions, one 
can find the homotopy deformation from the ODE limit cycle (procession) (25] [41] 
to the current PDF chaos as shown in Figure [T3] at the same parameter values. 

We also simulated the case of normal Gilbert damping a = 0.02. For all values 
of /32 € [0.01,0.3], the attractor is always non-chaotic. That is, the only attractor 
we can find is a spatially uniform limit cycle with small temporal oscillation as 
shown in Figure [TH 

Of course, when neither /3i nor (32 is zero, the bifurcation diagram is a combi- 
nation of the DC only and AC only diagrams. 
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Figure 13. Homotopy deformation of the attractors under differ- 



ent initial conditions. See http://www.math.missouri.edu/~'cli 



Figure 14. The attractor when a = 0.02, /32 = 0.21 and 
all other parameters' values are the same with Figure [121 See 
ihttp:/ /www. math. missouri.edu/~cli| 

6. Appendix: The Connection Between the Heisenberg Equation and 

the NLS Equation 

In this appendix, we will show the details on the connection between the ID 
cubic focusing nonlinear Schrodinger (NLS) equation and the Heisenberg equation 
p.ip . The nonhnear Schrodinger (NLS) equation 

(6.1) igt + fe + 2|gpg==0 , 
is a well-known integrable system with the Lax pair 

(6.2) 9,0 = (zA(73 + C/)0, 

(6.3) 5t0 = -(2aV3 + 2AJ7-l-T/)(/) , 
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where A is the complex spectral parameter, as is defined in ()2.5|) . and 

'^^{I'l)' "-I'l^-K",-. '5)- 

Lemma 6.1. Ifcj) — {(j)i,(j)2)'^ solves the Lax pair S6.2\) - ![K3\) at \, then {—(j)2, 4>i)'^ 
solves the Lax pair i6.S^) - [K^) at X. When q is even, i.e. q{^x) — q{x), then 
{(j>2{—x), (j)i{—x))'^ solves the Lax pair i6.S\} - [K^} at —A. When A is real, and cj) is 
a nonzero solution, then (— (/'2, is another linearly independent solution. For 

any two solutions of the Lax pair, their Wronskian is independent of x and t. 

For any real Ao, by the well-known Floquet theorem [28] and Lemma I6TT1 there 
are always two linearly independent Floquet (or Bloch) eigenfunctions (f>^ to the 
Lax pair (|6.2p - (|6.3p at A = Aq, such that 




{x + 27r) = p(t)+ (x) , {x + 27r) = p<j)- {x) , \pf = 1 ■ 

Since the Wronskian W{(f>'^ , <j)^) is independent of x and t, without loss of generality, 
we choose W{(t>'^ , cj)^) = 1. Then 

s=( 

is a unitary solution to the Lax pair at A = Aq: 

V -<P2 ^1 J 

Recall the definition of F ([27]) . let 

^ - ^ "^^-y -2^1^, \v2\'-\m\' ) ' 

i.e. 

mi + im2 ^ -2ipiip2 , ma = - |v?2p ■ 
Now for any solving the Lax pair (|6.2p - (|6.3p at A, define V' as 

Then ip solves the pair 

(6.4) V'. = *(A - Ao)FV; , 

(6.5) = -|2j(A2-A2)r+i(A-Ao)[r,r,]|^ . 
The compatibility condition of this pair leads to the equation 

(6.6) rt = -|4Aor, + ^[r,r,j| . 

Setting Ao = or performing the translation t = t, x = x — AX^t, equation (|6.6p 
reduces to the Heisenberg equation p.2p . Therefore, the Gauge transform S trans- 
forms NLS equation into the Heisenberg equation. Periodicity in x may not persist. 
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Example 6.2. Consider the temporally periodic solution of the NLS equation 

(EH), 

q ^ ae'^(*) , e{t) = 2ah + 7 . 
The corresponding Bloch eigcnfunction of the Lax pair (|6.2p - (|6.3p at A = is 

\e/2 



Then 







which is called a domain wall. 
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